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Theoretical Studies of Solar Lasers and Converters 
Iodine Laser Geometry. 

The geometry and setup for the n — CzF 2 I Iodine laser is illustrated in the Figure 1 . 
We are interested in the steady state laser operation with axial lasant flow. The following 
is a description of the mathematical modeling of this system. 

Chemical Kinetics. 

The chemical kinetics associated with this iodine laser has been disscused in the 
references 1 through 6 . Let 72 denote the radical n — C 3 F 7 the chemical kinetics can 
be summarized as follows: 

Photodissociation of RI and I 2 


0. 27 n m 

RI + hi/\ 72 4 - I* 

Cl 

0.5 pm 

I2 + hi/ 2 I + I* 

£2 


Laser action (stimulated emission) 


Quenching of I* (reaction rates have units [cm 3 /sec]) 

7 * + RAl + RI 

I*+h^I + h 

I* + 72-^7 + R 

I* + R 2 ^I + R 2 

I * + 7-^*7 + 7 


RI generation and two body recombinations (reaction rates have units [cm z / see]) 

72 + 7* -^+727 
72 + 7-^727 
72 + 72-^72 2 
72 + RI-^R 2 + 7 
72 + I 2 -^RI + I 
72 + RI^R 2 + 1* 

727 + 7*-^7 2 + 72 
7 + 727-^7 2 + 72 
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Ii formation (reaction rates have units [cm 6 /sec]) 


I* +1 + RI-^h + RI 
I + 1 + RI-^I 2 + RI 

r + i + i 2 -^i 2 + h 
i + i + h-^h + h 

I + 1 + r 2 -^i 2 + r 2 

I + I * + i ? 2 — “*-^2 + R2 


Wall reactions (reaction rates have units [cm 3 /sec]) 

Wall reactions have been proposed but have not as yet been implemented into the model. 
Two possible wall reactions are: 


I + I + wallah + wall 
R 2 -(■ / "h wall — — *R + RI + wall 


Values of parameters. 

The reaction rate coefficients in the above chemical reactions are assumed to have the 
bounds given in the Table 1. This Table also gives the final values used in the computer 
simulation. 


Differential equations for chemical kinetics. 

Using the notation [A] to denote the concentration of species A in units of cm -3 , the 
differential equations defining the chemical kinetics for the iodine laser can be expressed 
by the following set of coupled nonlinear differential equations involving the concentrations 
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[RI], [72], [R 2 ], [I 2 ], [I*] and [/]: 

& = kMl*} + k 2 [R][I } + k B [R][I 2 } - k 7 [r\[RI] - k 4 {R}[RI) 

- k 6 {R}[RI] - £ x [RI] + v 2 [R 2 \[I\ - k s [I}[RI] 

® - 6 [RI] - * 1 W1 - k 2 [R][I) - 2k 3 [R) 2 - k 4 [RI][R ] 

- k 6 [RI][R] - k s [R)\I 2 ] + v 2 [R 2 ][I] + k 7 [RI}[I*} + k 8 {I][RI } 

® = k 3 [R} 2 + k 6 \RI][R ] + k 4 [RI}[R) - v 2 [R 2 ][I } 

® = Cl [RI}[n[I] + c 2 [RI][I] 2 + c 3 [J 2 ] [/*][/] + c 4 [I 2 ][I} 2 

- £ 2 [7 2 ] + k 7 [RI)\r} - ks[R][I 2 \ + Vl [I} 2 
+ c 5 [I} 2 [R 2 \ + * 8 [*/][J] + C 6 [I}[I*}[R 2 } 

® = QyhlRI] + £ 2 [7 2 ] - k.lRWn - Q 2 [I 2 ]\r} 

- ccp([r] - i[J]) + ke[R}[RI] - $«[*][/*] - Q 4 [R 2 }[r\ 

- Qsini] - k 7 [RI][r } - c 6 [R 2 ][I*][I} 

® = uih] + Qiimn + QtMn - 2 c 5 [i} 2 [r 2 ] - m w 

+ cMin - \[I}) - cxiRini) -2 c 2 [RI][I) 2 - c 3 [I 2 ][r][I) 

- 2c 4 [7 2 ][7] 2 - k 2 [R}[I] + k 4 [Rl\[R] + Q»\r\[R] + <? 4 [/*p 2 ] 

+ Qs[r}[I} + * 5 [iZ][J 2 ] - V 2 [R 2 }[I] - 2v 1 [I} 2 - c 6 {R 2 ][r ][/] 

In the above differential equations we neglect compressibility and use the material deriva- 

<j \± H _ . ?[A }dz _ d[A .] d[A]^ 

dt dt dz dt dt dz 

where ^| = w is the flow rate in the axial direction. Note that the above system of nonlinear 
differential equations conserves the masses of the species involved in the reactions and at 
any point z we have the integrals 

[RI] + [J2] + 2[R 2 ] = constant 
[RI] + 2[I 2 ] + [/*] + [7] = constant 

In the above equations describing the reaction kinetics Q y is the quantum yield; £ 1 , £2 are 
the photodissociation rates which are dependent upon the solar simulator concentration 
co; o is the stimulated emission cross section; and p is the photon density in the optical 
cavity. As a first approximation we assume that 

6 = (3.04(10 -3 ))c o £ 2 = (3.38(10- 2 ))c 0 
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For the light flux density of lasing photons p we let p = p + + p- where p+ = p+{z, t ) 
denotes the photon density propagating in the positive z direction and p- — p- (z, t) de- 
notes the photon density propagating in the negative z direction. The differential equations 
for these photon densities are given by 


1 dp+ dp± 

c dt dz 
1 dp- dp- 

c dt dz 


= o' + (\r \ - j[/|) 

= <T / ,_([/*| -![/]) 


where c is the speed of light in the medium. In the above equations a[I*\ is the amplifi- 
cation factor resulting from population of the upper lasing level of the active medium and 
— ^<r[/]p+ is the decrease in photon density due to population of the lower lasing level. 

Initial and Boundary conditions. 

We consider the steady state solution to the above differential equations which are 
then subject to the initial condtions: 


[ RI ] = 3.5(10 16 )Po at z = 0 Pq is the pressure in torr 


with all of the other concentrations being zero at z = 0. 

Boundary conditions must be imposed upon the p + , p_ solutions due to the reflectiv- 
ities from the mirrors and Brewster windows of the geometry in Figure 2. In this geometry 
Ri, P 2 > w are given and Ruff,R 2 eff,t m are calculated. In this study Ruff = 0.95848 
and i? 2 e// = 0.81634 denote the effective reflectivities of the Brewster window-mirror com- 
bination where w = 0.98 denotes the transmission coefficient for the Brewster windows. 
The transmission coefficient of the output mirror is then given by t m = rv(l — R 2 ) where 
i ?2 is the reflectivity of the output mirror and is obtained from the relation 

Rieff = w 2 P 2 

This gives a transmission coefficient of t m = 0.147. 

At the faces z = 0 and z = L the differential equations for p+ and p- are subject to 
the boundary conditions: 


P+(°) = RleffP— (0) P-{L) = RieffP+{L) 

For the steady state , the solutions for p + and p- must satisfy the relation 

p+(,e)p_(z) = Kq = constant 
for all values of z. This condition requires that at z = 0 

0+(0)o-(0) = *..//[«-( 0)1 2 = -j^-MO)] 2 = Ko 

Ruff 
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and at z = L we must have 

„+(£)<>-(£) = = - r ’,//IM £)] 2 = *o 

These conditions require that 

P+( 0) = y/K^Tf and />+(£) = 

A shooting method is used to determine the initial conditions assigned to the steady state 
value of p+. For fixed values of Ruff and Rieff we quess at an initial value for Ko and 
integrate the system of steady state differential equations from z — 0 to z — L and compare 
the theoretical value of p+(L ) with the calcualted value. This iteration continues until this 
boundary condition is satisfied, in which case the power transmission is given by 

P = t»t m cp + (L) ( 1 W/cm 2 ) 



where e M is the energy quantum of lasing photons. 

Using the conversion factor that one solar constant equals 1.35 Kw/m 2 and given that 

12.9 Kw corresponds to 300 amps then the input power is assumed to be linear and 

T * r, \ 12.9 Kw 

Input power = I [amps) 

v ' 300 amps 


The output power is then given by 


output power = PdAf, 


where Pd is the power density and Ab is the cross sectional area of the laser beam. The 
power density varies with the beam diameter and the input power. For a beam diameter 
of d (cm) we have 


Input power 

Hd 


= Pd 


where t (cm) is the length of the tube section being radiated by the solar simulater. Using 
these relations we obtained good agreement between the theoretical and actual output 
power from the laser. 


Conclusion. The Figures 3 and 4 are obtained by holding Ruff constant and letting 
R^ef f vary. Thus, if the numerical integration stopped at some point z along the tube we 
determined the value of Ruf / necessary to satisfy the boundary conditions. For this value 
of i? 2e / / the transmission coefficient t m and resulting power density was calculated. 

Also a sensisitivity analysis was done on the parameters occuring in the differential 
equations describing the chemical kinetics. The results of this sensitivity study indicated 
that a changing of the coefficients 


ki, C 2 , Q i, A?2, Qz 

produced the largest changes in the power output. 
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Reaction 

coefficient 

Lower 

bound 

Nominal 

value 

Upper 
... bound 

Reference 

Value 

used 

Qi 

0.476E-16 

2.0E-16 

8.40E-16 

( 3 ) 

0.476E-16 

Q 2 

0.730E-11 

1.9E-11 

4.94E-11 

(3) (5) 

0.730E-11 

Q 3 

1.230E-18 

3.7E-18 

11.1E-18 

( 5 ) 

1.23E-18 

Q 4 

1.570E-16 

4.7E-16 

14.1E-16 

( 5 ) 

1.57E-16 

Q 5 

0.530E-14 

1.6E-14 

4.80E-14 

( 5 ) 

0.53E-14 

*1 

0.9E-13 

5.6E-13 

34.7E-13 

( 3 ) 

0.9E-13 

ki 

0.657E-11 

2.3E-11 

8.05E-11 

( 3 ) 

8.05E-11 

kz 

0.65E-12 

2.6E-12 

10.4E-12 

( 3 ) 

0.65E-12 

k 4 

1.0E-16 

3.0E-16 

9.0E-16 

(3) (5) 

1.0E-16 

k 5 

0.33E-11 

1.0E-11 

3.0E-11 

( 5 ) 

3.0E-11 

k 6 

1.0E-17 

3.2E-17 

10.24E-17 

( 3 ) 

1.0E-17 

k 7 

1.5E-19 

3.0E-19 

4.5E-19 

( 3 ) 

1.5E-19 

ka 

0.533E-23 

1.0E-23 

4.8E-23 

( 5 ) 

0.0 

Cl 

1.0E-33 

3.2E-33 

10.2E-33 

( 3 ) 

1.0E-33 

C2 

1.6E-32 

8.5E-32 

45.0E-32 

( 3 ) 

45.0E-32 

C3 

4.44E-32 

8.5E-32 

14.4E-32 

( 3 ) 

4.44E-32 

C 4 

2.92E-30 

3.8E-30 

4.94E-30 

( 3 ) 

4.94E-30 

C5 

3.6E-31 

4.8E-31 

6.0E-31 


3.6E-31 

Cq 

1.35E-32 

1.8E-32 

2.25E-32 


0.0 

Vi 

0.33E-12 

1.0E-12 

3.0E-12 


0.0 

v 2 

0.33E-11 

1.0E-11 

3.0E-11 


0.0 


Table 1. Parameter Values For Chemical Kinetics 
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FIGURE 1 

GEOMETRY FOR STEADY STATE LASER WITH AXIAL LASANT FLOW 




FIGURE 2 

MIRROR REFLECTIVITIES 
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THEORETICAL POWER OUTPUT (W / 
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APPENDIX 


COMPUTER PROGRAM FOR SIMULATION OF LASER POWER 
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PROGRAM NLAS7 ( I NPUT , OUTPUT , TAPE5= I NPUT , TAF'E6=0UTPUT , T APE8 ) 

MAIN PROGRAM 

THIS VERSION CONTAINS AN AUTOMATICALLY CHOSEN VARIABLE STEP 
SIZE. IT ALSO CONTAINS WALL EFFECT REACTIONS. 

COMMON /BLK2/ K1 , K2, K3, K4, K5, K6, K7, K8, Cl , C2, C3, C4, 01 , Q2, A, AD, VI , V2 
1 , GG 

COMMON /BLK3/ B , B2 , B3 , C , AOO , BOO , EPSNU , OMEGA , C5 , C6 , Q3 , Q4 , Q5 

COMMON /BLK4/ CHSI 10, CHSI20, ABARO, Z 1BAR, LC 

COMMON /BLK7/ ABC, COO, CO, 0MEG1 , P, R1 , R2, TM 

COMMON /BLK8/ ZOL 

REAL LC , K 1 , K2 , K3 , K4 , K5 , K6 , K7 , K8 

WRITE (6,10) 

10 FORMAT (IX, ’ START OF PROGRAM ’) 

DEFAULT VALUES 

ZOL = 5.70 
P = 40. 

LC = 6.0 
COO = . 5E25 
R1 = .9975 
R2 = .975 
A = 0.0 
TM = 1-R2 
0MEG1 = 5000. 

CON = 2.0E4 

NAMELIST /F'ARAM/ P, 0MEG1 , CON, COO, R1 , R2, TM, LC, ZOL, A 
20 CONTINUE 

READ (5, PARAM) 

IF (EOF (5)) GO TO 30 
GO TO 50 
30 WRITE (6,40) 

40 FORMAT (IX, ’ END OF FILE ENCOUNTERED-STOP’ ) 

STOP 1313 


P=PRESSURE, TORR 
0MEG1=FL0W RATE, CM/SEC 
CON=PEAK CONCENTRATION 
C00= INITIAL GUESS AT RHO-PLUS AT ZERO 
WHICH IS SQUARE OF (C00*R1> 

Rl= REFLECTIVITY AT LEFT END 

R2= REFLECTIVITY AT RIGHT END 

TM= TRANSMISSION COEFFICIENT (OUTPUT MIRROR) 

ZOL=POINT ALONG AXIS WHERE MAXIMUM ILLUMINATION OCCURS 
IN THE CASE ILLUMINATION IS A BELL SHAPED CURVE 
IN THE CASE OF A SQUARE WAVE, 2*Z0L IS CUT OFF POINT 
THE POINT 2*Z0L IS WHERE ILLUMINATION CUTS OFF 
LC=LENGTH OF CAVITY 
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50 CONTINUE 

CM IN = 1.0E10 
CMAX = 1.0E30 
CO = CON 
CU = CON 
OMEGA = 0MEG1 
WRITE (6,60) 

60 FORMAT (///) 

WRITE (6,70) ZOL, CON, OMEGA, COO, R1,R2,P 
70 FORMAT <1X,T30,’Z0L - ’ , E15. 7, T55, ’ CON = ’ , E15. 7, T80, ’ OMEGA = ’, 

1 E15.7, /, IX, T5, 'COO = ’ , E15. 7, T30, ’ R1 = ’ , F10. 7, T55, ’ R2 = 

2 F10. 7, T81, ’ P =’ , E15. 7) 

SET UP COEFFICIENTS IN DIFFERENTIAL EQUATIONS 
CALL COEFFS 

INTEGRATE DIFFERENTIAL EQUATIONS FROM Z=0 TO Z=LC 

XI = COO 
CALL INTEG 

INTERVAL HALVING SCHEME 

W1 AND W2 ARE WEIGHTS FOR INTERVAL HALVING SCHEME FOR 
DETERMINING COO WHICH SATISFIES BOUNDARY CONDITIONS 

Y1 = ABC 

IF (Yl.LT.O) PER = . 1 
IF(Yl.LT.O) PER=. 9 
IF (Yl.GT.O) PER = 10. 

IF(Yl.GT.O) PER=1 . 1 
80 CONTINUE 

EITHER INCREASE OR DECREASE INITIAL GUESS TO COO 
COO = (PER) *C00 
STOP IF COO EXCEEDS BOUNDS 
IF (COO.LT.CMIN) STOP 5432 
IF (COO. GT. CMAX) STOP 2345 
X2 = COO 
; CALL INTEG 
Y2 = ABC 

IF ( (Yl*Y2).LT.O) GO TO 90 
XI = COO 
Y 1 = Y2 
GO TO 80 
90 CONTINUE 

W1 AND W2 ARE WEIGHTING CONSTANTS FOR INTERVAL HALVING 


Y 1 AND Y2 ARE OF OPPOSITE SIGN 
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W3 = ABS (Y1 ) +ABS (Y2) 

W1 = ABS<Y1)/W3 
W2 = ABS(Y2)/W3 
COO = W2*X1+W1*X2 
100 CONTINUE 
CALL INTEG 
X3 = COO 
Y3 = ABC 

IF (ABS(Y3) .LT. 0.001) GO TO 20 
IF < <Y1*Y3) .LT.O) GO TO 110 

Y1 & Y3 ARE OF THE SAME SIGN 

XI = X3 
Y 1 = Y3 

W3 = ABS < Y1 ) +ABS < Y2) 

W1 = ABS(Y1)/W3 
W2 = ABS(Y2)/W3 
COO = W2*X 1+W1*X2 
GO TO 100 
110 CONTINUE 

Y 1 & Y3 ARE OF OPPOSITE SIGN 


X2 = X3 
Y2 = Y3 

W3 = ABS (Y1 ) +ABS (Y2) 
W1 = ABS(Y1)/W3 
W2 = ABS(Y2)/W3 
COO = W2*X1+W1*X2 
GO TO 100 

STOP 1234 
END 


FUNCTION CHS 1 1 (Z) 


COMMON /BLK4/ CHSI 10, CHSI20, ABARO, Z 1BAR, LC 
COMMON /BLK8/ ZOL 

REAL K 1 , K2 , K3 , K4 , K5 , K6 , K7 , KB , Q3 , Q4 , Q5 , LC 
IF (Z.LT. ABARO) GO TO 10 
IF (Z.LT.Z1BAR) GO TO 20 

Z GREATER THAN Z1BAR 

10 CHSI 1 = 0.0 
RETURN 
20 CONTINUE 

CHSI 1 = CHSI 10 

RETURN 

END 
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FUNCTION CHS I 2 <Z> 


COMMON /BLK4/ CHSI 10, CHSI20, ABARO, Z 1BAR, LC 
COMMON /BLK8/ ZOL 

REAL K 1 , K2 , K3 , K4 , K5 , K6 , K7 , K8 , Q3, Q4, Q5, LC 
IF (Z.LT. ABARO) GO TO 10 
IF (Z.LT.Z1BAR) GO TO 20 

Z GREATER THAN Z1BAR 

10 CHS I 2 = 0.0 
RETURN 
20 CONTINUE 

CHS I 2 = CHS I 20 

RETURN 

END 


SUBROUTINE COEFFS 

THIS SUBROUTINE DEFINES THE COEFFICIENTS IN THE DIFFERENTIAL 
EQUATIONS TO BE SOLVED. 

IMPLICIT REAL*8 ( A-H, K, L, 0-Z) 

COMMON /BLK2/ K1 , K2, K3, K4 , K5, K6, K7, K8, Cl , C2, C3, C4, Q1 , Q2, A, AD, VI , V2 
1 , GG 

COMMON /BLK3/ B, B2, B3, C, AOO, BOO, EPSNU, OMEGA, C5, C6, Q3, Q4, Q5 
COMMON /BLK4/ CHSI 10, CHSI20, ABARO, Z 1BAR, LC 
COMMON /BLK7/ ABC, COO, CO, 0MEG1 , P, R1 , R2, TM 
COMMON /BLK8/ ZOL 

REAL K1,K2,K3,K4,K5,K6,K7,K8,Q3,Q4,Q5,LC 

COEFFICIENTS IN THE DIFFERENTIAL EQUATIONS 

OMEGA = 0MEG1 
ABARO =0.0 

ABARO= START OF ILLUMINATION 
Z1BAR=2*Z0L = POINT ON AXIS WHERE ILLUMINATION STOPS 

CHSI 10 = (3. 04E-3) *C0 
CHS I 20 = <3. 38E-2) *C0 
Z1BAR = 2*Z0L 
EPSNU = 1.5E-19 

WATTS*SEC 

AOO = 2.0E17 
BOO = .443 
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ccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

C (CHS 1 1 ) HV + RI — > R + I* 

C (CHSI2) HV + 12 — > I* + I 

C 


c 

(Kl) 

I* + 

R — > RI 

c 

(K2) 

I + 

R — > RI 

c 

(K3> 

R + 

R --> R2 

c 

( K4 ) 

R + 

RI — > R2 + I 

c 

(K5> 

R + 

12 — > RI + I 

c 

<K6> 

R + 

RI — > R2 + I* 

c 

(K7> 

I* + 

RI — > 12 + R 

c 

(K8> 

I + 

RI — > 12 + R 

c 




c 

(Ql) 

I* + 

RI — > I + RI 

c 

(Q2) 

I* + 

12 — > I + 12 

c 

(Q3) 

I* + 

R — > I + R 

c 

<Q4) 

I* + 

R2 — > I + R2 

c 

(Q5) 

I* + 

I --> I + I 

c 




c 

(Cl) 

I* + 

I + RI — > 12 

c 

(C2) 

I + 

I + RI — > 12 

c 

(C3) 

I* + 

I + 12 — > 12 

c 

(C4) 

I + 

I + 12 — > 12 

c 

(C5) 

I + 

I + R2 — > 12 

c 

(C6) 

I + 

I* + R2 --> 12 


i- Rl 
RI 

i- 12 
12 
R2 
i- R2 


WALL REACTIONS 


C 
C 
C 
C 

C (VI) I + I + WALL — > 

C ( V2) R2 + I + WALL — > 

C 

ccccccccccccccccccccccccccccccccccccccccccccccccc 

c 


12 + WALL 
R + RI + WALL 


C 

Kl 

as 

34.7E-13 

C 

Kl 

SB 

5.6E-13 

C 

Kl 

= 

. 9E— 13 

C 

K2 

=r 

S.05E-11 

C 




C 

K2 

= 

2.3E-11 

C 

K2 

= 

.657E-11 

C 




C 

- K3 

3S 

10.4E-12 

C 

K3 

=r 

2.6E-12 

C 

K3 

_ 

. 65E-12 

C 




C 

K4 

= 

9.E-1 

C 

K4 

— 

3.E-16 

C 

K4 

_ 

l.E-16 



non o noon nnno noon onoo nnno nnno nonnn nnnn oononon 


K.5 


3.E-11 


K5 

= 

1.0E-11 

K5 

= 

.33E-11 

K6 

sz 

10.24E-17 

K6 

= 

3.2E-17 

K6 

3= 

1 . E— 17 

K7 

sr 

4.5E-19 

K.7 

= 

3.0E-19 

K7 

s 

1.5E-19 

K8 

32 

4.8E-23 

K8 

= 

1 . 6E-23 

K8 

= 

. 533E-23 

K8 

= 

0.0 

Q1 

ss 

8.4E-16 

Q1 

= 

2.0E-16 

Q 1 

s 

. 476E-16 

Q2 

r : 

4.94E-U 

Q2 

= 

1.9E-11 

Q2 

= 

. 73E-1 1 

Q3 

sr 

11. IE-18 

Q3 

= 

3.7E-18 

03 

= 

1.23E-18 

Q4 

32 

14. IE-16 

04 

= 

4.7E-16 

Q4 

= 

1.57E-16 

0.5 

32 

4.8E-14 

05 

32 

1.6E-14 

05 

= 

. 53E-14 

Cl 

=2 

10.2E-33 

Cl 


3.2E-33 

Cl 

= 

1 . E— 33 

C2 

= 

45.E-32 

C2 

32 

8.5E-32 

C2 

23 

1 . 6E-32 
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C3 


14. 4E 

-32 

C3 

= 

8.E-32 

C3 

= 

4.44E 

CM 

w 

1 

C4 

= 

4.94E 

-30 

C4 

= 

3.8E- 

30 

C4 

= 

2.92E 

-30 

C5 

— 

6.E-3 

;i 

C5 


4.8E- 

31 

C5 

= 

3.6E- 

31 

C6 

- 


2.25E-32 

C6 



1.8E-32 

C6 



1 . 35E-32 

C6 

= 

o 

• 

o 


VI 

- 


3.0E-12 

VI 

= 


1.0E-12 

VI 

= 


. 33E- 1 2 

VI 

= 

o 

■ 

o 


V2 



3.0E-U 

V2 

= 


1.0E-11 

V2 

= 


.33E-11 

V2 

= 

o 

• 

o 



Kl = . 903E-13 
K2=80. E-12 
K3=. 65E-12 
K4=l . 000E-16 
K5=3 . 009E- 1 1 
K6=l . OE-17 
K7=. 1517E-18 
K8= 1 . 6E-23 


Ql = . 476E-16 
Q2=1.9E-11 
Q3=. 1235E-17 
Q4=l . 57E-16 
05= . 53E-14 


Cl=l . 053E-33 
C2=45. OE-32 
C3=. 4447E-31 
C4=4. 94E-30 
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non 


C C5=3. 6E-31 

C C6= 1 . 8E-32 

C 

C Vl = 1.0E-12 

C V2= 1.0E-11 

C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

C 

GG = 2 * ( . 18/LC) **2 
C = 3.0E10 
B = P* (3. 5E16) 

B2 = B*B 
B3 = B2*B 

WRITE OUT COEFFICIENTS 

WRITE (6,10) K1,K7,Q1,C5 
WRITE (6,20) K2, K8, Q2, VI 
WRITE (6,30) K3,C1,Q3,V2 
WRITE (6,40) K4,C2,Q4 
WRITE (6,50) K5,C3,Q5 
WRITE (6,60) K6,C4,C6 

10 FORMAT (T5, ’ K1 = ’ , E15. 7, T30, ’ K7 = ’ , E15. 7, T60, ’ Q1 = ',E15.7,TB5, 
1 ' C5 = ’ , E15. 7) 

20 FORMAT (T5,’K2 = ’ , E15. 7, T30, ’ K8 = ’ , E15. 7, T60, ’ Q2 = ',E15.7,TB5, 
1 'VI = ’ , E15. 7) 

30 FORMAT (T5,'K3 = ’ , E15. 7, T30, ’ Cl = ' , E15. 7, T60, ' Q3 = ',E15.7,T85, 
1 ' V2 = ’ ,E15.7) 

40 FORMAT (T5,’K4 = ’ , E15. 7, T30, ’ C2 = ' , E15. 7, T60, ’ Q4 = ',E15.7) 

50 FORMAT (T5,'K5 = ' , E15. 7, T30, ’ C3= , ’ , E15. 7, T60, ’ Q5 = ',E15.7) 

60 FORMAT <T5,’K6 = ’ , E15. 7, T30, ' C4= , ’ , E15. 7, T60, ' C6 = ’,E15.7) 
RETURN 
END 
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SUBROUTINE FUN <N,Z,Y,F> 

THIS SUBROUTINE DEFINES THE RISHT HAND SIDE 
OF THE DIFFERENTIAL EQUATIONS FOR THE CHEMICAL KINETICS 
IMPLICIT REAL*8 ( A-H, K, L, O-Z) 

DIMENSION Y<7) ,F(7) 

COMMON / BLK 1 / X 7, POWER 
EXTERNAL CHSI1,CHSI2 

COMMON /BLK2/ K 1 , K2 , K3 , K4 , K5 , K6 , K7 , K8 , C 1 , C2 , C3 , C4 , Q 1 , Q2 , A , AD , V 1 , V2 
1 ,68 

COMMON /BLK3/ B, B2, B3, C, AOO, BOO, EPSNU, OMEGA, C5, C6, Q3, Q4, Q5 

COMMON /BLK4/ CHS 1 10, CHS 1 20, ABARO, Z 1BAR, LC 

COMMON /BLK7/ ABC, COO, CO, 0MEG1 , P, R1 , R2, TM 

REAL K 1 , K2 , K3 , K4 , K5 , K6 , K7 , K8 , Q3 , Q4 , Q5 , LC 

QY=QUANTUM YIELD 

QY = 1.0 

F< I), 1=1,6 ARE RATES OF CHANGES FOR THE CONCENTRATIONS 
F<1)=DCRI]/DZ F<2)=D£R3/DZ 

F <3) =DCR23/DZ F (4) =DC 123 /DZ 

F<5)=DCI*]/DZ F (6) =DC I ] /DZ 

X8 = C00/(Y(7)*B) 

X7STAR = Y (7) *B+X8 
DIF = Y (5) - . 5*Y (6) 

CALL SIGMA (SIG2) 

SIG = SIG2 

F ( 1 ) = K1*B*Y (2) *Y (5) +K2*B*Y <2) *Y (6) -CHSI 1 (Z) *Y < 1 ) -K4*B*Y < 1 ) *Y (2) + 

1 K5*B*Y ( 2 > * Y ( 4 ) -K7*B*Y < 5 ) *Y < 1 ) -K6*B*Y ( 2 ) * Y ( 1 ) +V2*B*Y ( 3 ) * Y ( 6 ) -K8* 

2 B*Y<6)*Y<1) 

F (2) = CHSI 1 ( Z) *Y ( 1 ) -K1*B*Y (2) *Y <5) -K2*B*Y (2) *Y<6> -2*K3*B*Y (2) *Y <2 

1 ) -K4*B*Y < 1 ) *Y (2) -K6*B*Y ( 1 ) *Y (2) -K5*B*Y (2) *Y <4) +V2*B*Y (3) *Y <&) + 

2 K7*B*Y<5)*Y<1)+KB*B*Y<6>*Y<1> 

F (3) = K3*B*Y (2) *Y (2) +K6*B*Y ( 1 ) *Y <2) +K4*B*Y ( 1 ) *Y (2) -V2*B*Y <3) *Y (6) 
A1 = C1*B2*Y<1)*Y(5)*Y(6)+C2*B2*Y<1>*Y<6>*Y<6>+ 

1 C3*B2*Y<4)*Y(5)*Y<6) 

A2 = C4*B2*Y (4) *Y (6) *Y (6) -CHS 1 2 ( Z ) *Y (4) +K7*B*Y <5) *Y ( 1 ) - 
1 K5*B*Y (2) *Y (4) +V1*B*Y ( 6 ) *Y <6) +C5*B2*Y (6) *Y <6> *Y (3) 

F<4> = A1+A2+K8*B*Y (6) *Y ( 1 ) +C6*B2*Y (6) *Y <5) *Y (3) 

A3 = QY*CHSI1 (Z)*Y(1)+CHSI2<Z)*Y(4)-K1*B*Y(2)*Y<5> 

A4 = -C1*B2*Y ( 1 ) *Y (5) *Y (6) -C3*B2*Y (4) *Y (5) *Y (6) -Q1*B*Y ( 1 ) *Y (5) 

A5 = -Q2*B*Y(4)*Y(5)-C*SIG*X7STAR*DIF+K6*B*Y<2)*Y(1> 

F (5) = A3+A4+A5-Q3*B*Y (5) *Y (2) -Q4*B*Y (5) *Y <3> -Q5*B*Y (5) *Y (6) -K7*B* 
1 Y(5)*Y<1)-C6*B2*Y(6>*Y<5)*Y<3> 

A6 = CHS 12 ( Z ) *Y (4) +Q1*B*Y < 1 ) *Y (5) +Q2*B*Y <4> *Y <5) -2*C5*B2*Y (6) *Y (6) 

1 *Y (3) -K8*B*Y <6) *Y ( 1 ) 

A7 = C*SIG*X7STAR*DIF-C1*B2*Y(1)*Y<5)*Y<6> 

A8 = -2*C2*B2*Y ( 1 ) *Y (6) *Y (6) -C3*B2*Y (4) *Y <5) *Y (6) 

A9 = -2*C4*B2* Y < 4 ) *Y < 6 ) * Y ( 6 ) -K2*B*Y ( 2 ) *Y ( 6 ) +K4*B*Y ( 1 ) *Y ( 2 ) 

A10 = Q3*B*Y (5) *Y <2) +Q.4*B*Y (5) *Y (3) +Q5*B*Y (5) *Y (6) +K5*B*Y (2) *Y (4) - 
1 V2*B*Y <3> *Y (6) -2*V1*B*Y (6) *Y (6) 

F<6) = A6+A7+A8+A9+A10-C6*B2*Y<6)*Y<5)*Y(3) 

DO 10 I » 1, 6 

F < I ) = F (I) /OMEGA 
10 CONTINUE 

F (7) = Y (7) *DIF*B*SIG 

RETURN 

END 
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SUBROUTINE SIGMA (SIG) 

THIS SUBROUTINE DEFINES THE CROSS SECTION SIGMA 

COMMON /BLK3/ B, B2, B3, C, AOO, BOO, EPSNU, OMEGA, C5, C6, Q3, Q4, Q5 

COMMON /BLK7/ ABC, COO, CO, OMEG1 , P, R1 , R2, TM 

REAL NU , NUS , NUO , NU 1 , NU2 , NU3 , NU4, NU5 

PI = 3. 14159 

PIS = PI*PI 

NU = C/1.315246E-4 

NUS = NU*NU 

PISNUS = PIS*NUS*4. 

6 = 0 
CS = C*C 
NUO = NU 

NU1 = NUO+. 141*C 

NU2 = NU1+.068*C 

NU3 = NUO-. 427*C 

NU4 = NU3- . 026*C 

NUS = NU4-.068*C 

DELTA23 = NU-NU5 

DELTA22 = NU-NU4 

DELTA21 = NU-NU3 

DELTA34 = NU-NUO 

DELTA33 = NU-NU1 

DELTA32 = NU-NU2 

TEMPO = 293 

TWALL = TEMPO 

T 1 = TWALL 

A = 5.434 

A1 = A*2. 4/7. 7*CS 

A2 = A*3. 0/7. 7*CS 

A3 = A*2. 3/7 . 7*CS 

A4 = A*5. 0/7 . 7*CS 

A5 = A*2.2/7.7*CS 

Ad> = A*0.6/7.7*CS 

FUGTEMP * SORT (Tl/300. ) 

ALPHAM = 1 . 88E7*FUGTEMP 
DELDOP = 2.51 E8*FUGTEMP 
DELNU = DELDOP+ALPHAM*P 

SIGMA23 = Al/ <PISNUS*DELNU) / ( 1+ (2. *DELTA23/ DELNU) **2) *5. /12 
SIGMA22 = A2/<PISNUS*DELNU)/<1+(2.*DELTA22/DELNU>**2>*5./12 
SIGMA21 = A3/ (PISNUS*DELNU) / < 1+ (2. *DELTA21 /DELNU) ** 2 ) *5. /12 
SIGMA34 = A4/(PISNUS*DELNU)/<1+<2.*DELTA34/DELNU>**2>*7./12 
SIGMA33 = A5/ (PISNUS*DELNU) / ( 1+ (2. *DELTA33/DELNU) **2) *7. / 12 
SIGMA32 = A6/ (PISNUS*DELNU) / ( 1+ (2. *DELTA32/ DELNU) **2) *7. / 12 
SIGMAT = S I GMA23+S I GMA22+S I GMA2 1 +S I GMA34+S I GMA33+S I GM A32 
SIG = SIGMAT 
RETURN 
END 
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SUBROUTINE INTEG 


C 

C THIS SUBROUTINE INTEGRATES THE SYSTEM OF DIFFERENTIAL EQUATIONS 

C USING A VARIABLE STEP SIZE 7TH ORDER RUNGE KUTTA METHOD. 

C IMPLICIT REAL*8<A-H,K,L,0-Z) 

C 

DIMENSION Y<7) ,F<7) , Y0<7) , X (7) , WK<49> 

COMMON /BLK1 / X 7, POWER 

COMMON /BLK2/ K1 , K2, K3, K4, K5, K6, K7, K8, Cl , C2, C3, C4, Q1 , Q2, A, AD, VI , V2 
1 ,GG 

COMMON /BLK3/ B, B2, B3, C, AGO, BOO, EPSNU, OMEGA, C5, C&, Q3, Q4, Q5 
COMMON /BLK4/ CHS 1 10, CHS 1 20, ABARO, Z 1BAR, LC 
COMMON /BLK7/ ABC, COO, CO, 0MEG1 , P. R1 , R2, TM 
EXTERNAL FUN, CHS 1 1 , CHS 1 2 

REAL K 1 , K2 , K3 , K4 , K5 , K6 , K7 , K8 , Q3 , Q4 , Q5 , LC 
C 

C INTEGRATE SYSTEM FROM Z=0 TO Z=LC USING RUNGE-KUTTA METHOD 

C 

C X < 1 ) =RI 

C X (2) =R 

C X (3) =R2 

C X (4) =12 

C X (5) =1* 

C X (6) =1 

C X (7) =RHO+ 

C X8=RH0- 

C X9=I*-.5*I 

C 

WW2 = < . 98) **2 
N = 7 

TOL = 1.0E-6 
PD = 1.0 
MTH = 1 
H = 1.00 

HMIN = H/ 10000000. 

HMAX = H/ 100. 

HUSE = HMIN* 10 
I ERR = 0 
C 

C INITIAL CONDITIONS 

C 

ZO = 0.0 
YO < 1 ) = 1.0 
Z1 = 0.0 
DO 10 I = 2, 6 
YO ( I ) =0.0 
10' CONTINUE 
C 

C GUESS AT INITIAL CONDITIONS FOR X<7) AND X<8) 

C INITIAL CONDITION FOR X70=RH0+ DEPENDS UPON COO VALUE 

C 

X70 = SQRT (R1*C00) 

YO (7) = X70/B 
C 
C 
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C PRINT HEADER 

WRITE (6,20) 

20 FORMAT <///,T7, ’ Z’ ,T20, ’ CRI3’,T32, ’ CR3, ’ ,T45, ’ CR23, ’ ,T57, ’ CI23 

1 T69, ’ CI*r ,TBO, ’ [I] ’,T91,'CRH0+3’,T103,’CRH0-3’ ,T112, 

2 ’ INVERSION" ) 

30 CONTINUE 

C 

c 

DO 40 I = 1, 7 
X(I) = B*YQ ( I ) 

40 CONTINUE 
C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

C 

C USE CONSERVATION LAWS THE ABOVE SYSTEM OF DIFFERENTIAL 

C EQUATIONS HAS THE TWO INTEGRALS 

C 

C CRI 3 + CR] +2CR23 = CONSTANT= B 

C OR 

C Y<1) + Y<2) + 2Y(3)= 1.0 

C 

C AND CRI 3 + 2CI23 + CI*3 + CI3 = B 

C OR 

C Yd) + 2Y (4) + Y (5) + Y (6) =1.0 

C 

Y0TH2 = 1 . O-YO (1 ) -2. *Y0 (3) 

Y0TH6 = 1 . O-YO ( 5 ) -2 . 0*Y0 ( 4 ) -YO ( 1 ) 

C CHECK. ON INTEGRATION ROUTINE COMPARED TO ABOVE INTEGRALS 

C TEST TO SEE IF Y0TH2 IS ZERO 

C 

ERROR 1 =0.0 

IF (Y0TH2.EQ.0.0) GO TO 50 
ERROR 1 = (YO (2) -Y0TH2) *100/Y0TH2 
50 CONTINUE 

ERR0R2 = 0.0 
C 

C TEST TO SEE IF Y0TH6 IS ZERO 

C 

IF (Y0TH6.EQ.0.0) GO TO 60 
ERR0R2 = (YO (6) -Y0TH6) *100/Y0TH6 
60 CONTINUE 

IF ( ( ERROR 1. OR. ERR0R2) .GT. .5) WRITE (6,70) 

70 FORMAT (IX, ’ ERROR 1 OR ERR0R2 IS OUT OF BOUNDS ’) 

C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

C END OF TEST FOR CONSERVATION LAWS BEING SATISFIED 
C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

X8 = COO/X (7) 

X9 = X (5) -. 5*X (6) 

X7STAR = X (7) +X8 
C 
C 
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C USE SUBROUTINE SISMA TO CALCULATE CROSS SECTION SIGMA 

CALL SIGMA <SIG2> 

PRINT OUT VALUES FOR CONCENTRATIONS , RHO+, RHO- AND INVERSION 

WRITE (6,80) ZO, (X (I) , 1=1,7) , X8, X9 
80 FORMAT ( 1 X , E12. 5, 8E12. 5, E12. 5, E12. 5) 

USE 7TH ORDER RUNGE KUTTA INTEGRATION SCHEME WITH VARIABLE STEP 
STEP SIZE CAN VARY FROM HMIN TO HMAX 
ZO IS STARTING VALUE FOR Z 

Z1 IS NEXT STOPPING POINT IN INTEGRATION SCHEME 
TOL IS TOLERENCE 

I ERR IS ERROR CODE TO DETERMINE IF INTEGRATION WAS SUCCESSFUL 
Z1 = Zl+H 

CALL RKF7 < N , ZO , Z 1 , YO , TOL , FUN , PD, MTH, HMIN, HMA X , HUSE , WK , I ERR ) 

IF (IERR.NE.O) WRITE (6,90) 

IF (IERR.NE.O) STOP 1717 

90 FORMAT (IX,’ I ERR IS NOT ZERO ’) 

X (7) = B*YQ (7) 

X8 - COO / X ( 7 ) 

IF ( (Z1+.5*H> .GE. (LC+.2) ) GOTO 110 

GO TO 30 
DO 100 I = 1, 7 
X(I) = B*YO ( I ) 

100 CONTINUE 

X8 = COO/X (7) 

X9 = X (5) - . 5*X (6) 

110 CONTINUE 

XX7L = B*YO (7) 

RHOPL = X70/ (SORT (R2*R1 ) ) 

RHOPL=RHO-PLUS AT Z=L THEORETICAL VALUE 
XX7L= CALCULATED VALUE OF RHO-PLUS AT Z= 

ABC =D I FFERENCE=X X7L-RH0PL 

DIF = (((XX 7L-RH0PL ) / RHOPL ) * 1 00 ) 

ABC = DIF 

WRITE (6,120) DIF, RHOPL, XX7L, COO 
120 FORMAT ( 1 X , ' DIFFERENCE = ’ , E18. 9, 2X, ’ THEORETICAL=’ , E18. 9, 2X , 

1 ' ACTUAL « ’ , E18. 9, 2X, ’ COO = ’,E18.8) 

XX7L = X (7) 

TM = 1-R2/ (WW2) 

IF (TM.LT.O) GO TO 30 
- POWER = EPSNU*TM*C*XX7L 
WRITE (6,130) R1,R2, POWER, TM,Z1,P 

130 FORMAT (1X,’R1 = ' , F10. 7, 1 X, ’ R2 = ’ , F10. 7, IX, ’ POWER =’ , E18. 10, 1 X , 
1 ’ TM = ’ ,F10.8, IX, ’L = ’^15.7,’? = '.F15.7) 

IF ( (Z1+.5*H) .GE. (LC+.2) ) GOTO 140 
GO TO 30 
140 CONTINUE 
RETURN 
END 
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